Os índigenas Pima - diabetes

Gabriel de Jesus Pereira

2025-04-05

Sobre o banco de dados

Sobre o banco de dados

  • Foi originado do National Institute of Diabetes and Digestive and Kidney Diseases.

  • Essa base foi criada com base nos Pima, um grupo de nativos americanos que vivem em uma área que atualmente abrange o centro e o sul do estado do Arizona.

  • A base possui 768 observações e 9 colunas. O objetivo é classificar se um paciente tem diabetes com base em algumas características do paciente.

Descrição das colunas

  • Pregnancies: Número de gestações da paciente.

  • Glucose: Concentração de glicose plasmática em teste oral de tolerância à glicose.

  • BloodPressure: Pressão arterial diastólica (mm Hg).

  • SkinThickness: Espessura da dobra cutânea do tríceps (mm).

  • Insulin: Nível sérico de insulina (mu U/ml).

  • BMI: Índice de Massa Corporal (peso em kg / altura em m², IMC).

  • DiabetesPedigreeFunction: Histórico familiar de diabetes (mede a predisposição genética).

  • Age: Idade do paciente.

  • Outcome: 0 é não diabético e 1 é diabético.

Análise exploratória de dados

Quantidade de valores ausentes

library(ranger)
library(dplyr)
library(ggplot2)
library(kknn)
library(naivebayes)
library(glmnet)
library(discrim)
library(tidyverse)
library(purrr)
library(tidymodels)
library(visdat)
library(kernlab)
library(corrplot)
library(baguette)

# data <- read_csv("aprendizagem_maquina/diabetes.csv") |> mutate(Outcome = as.factor(Outcome))
data <- read_csv("diabetes.csv") |> mutate(Outcome = as.factor(Outcome))
cols <- c("Glucose", "BloodPressure", "SkinThickness", "Insulin", "BMI")
data[cols][data[cols] == 0] <- NA


Ausentes <- data |>
  vis_dat() +
  labs(
    y = "Observações"
    ) +
  scale_fill_manual(
    values = viridis::mako(10)[5:6],
    labels = c("Qualitativa", "Quantitativa", "Ausentes")
    ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
  )

Quantidade de valores ausentes

Apenas as colunas de glicose, pressão arterial, espessura da dobra cutânea do tríceps, nível de insulina e IMC. Aquela que mais possui valores ausentes é a variável de nível de insulina.

Frequência

set.seed(123)
data_split <- initial_split(data, prop = 0.8, strata = Outcome)
train_data <- training(data_split)
test_data <- testing(data_split)

rec <- recipe(Outcome ~ ., data = train_data) %>%
  step_impute_median(all_numeric_predictors())

prep_rec <- prep(rec)
train_data_eda <- bake(prep_rec, new_data = train_data)
test_data_eda <- bake(prep_rec, new_data = test_data)

Frequência <- train_data_eda |>
  ggplot(aes(x = forcats::fct_infreq(as.factor(Outcome)))) +
  geom_bar(
        aes(
          y = after_stat(count) / sum(after_stat(count)),
          fill = Outcome
          ),
        show.legend = FALSE,
  ) +
  theme_bw() +
  theme(axis.title.x = element_blank()) +
  scale_fill_manual(
    values = viridis::mako(5)[3:5]
    ) +
  labs(y = NULL, x = NULL)

Frequência

A base de dados possui um problema de desbalanceamento das classes. A classe com maior frequência é de não diabéticos, com pouco mais de 60%.

Distribuição

Distribuição <- train_data_eda |>
  pivot_longer(
    !Outcome,
    names_to = "variable",
    values_to = "values"
  ) |>
  mutate(Outcome = as.factor(Outcome)) |>
  ggplot(
    aes(
      x = values,
      fill = Outcome
      )
    ) +
  geom_violin(
    aes(y = variable),
    draw_quantiles = c(0.25, 0.5, 0.75),
    ) +
  labs(x = "", y = "") +
  facet_wrap(vars(variable), scales = "free") +
  scale_fill_manual(
    values = viridis::mako(5)[3:5],
    name = "Outcome"
    ) +
  theme_bw() +
  theme(axis.text.y = element_blank())

Distribuição

As variáveis apresentam uma assimetria negativa, sendo a variável de pressão arterial a única mais simétrica. Além disso, observando a distribuição dos diabéticos e não diabéticos para cada variável, eles apresentam um padrão de comportamento bastante similar.

Correlação

data_cor <- train_data_eda |>
  (\(x) {
     scaled <- x |>
       select_if(is.numeric) |>
       scale() |>
       as_tibble()

     x |>
       select_if(is.character) |>
       cbind(scaled)
  })()

data_cor |>
  select_if(is.numeric) |>
  cor() |>
  corrplot(
      method="circle",
      col = viridis::mako(200),
      diag = FALSE,
      addgrid.col = "#00000040",
      tl.pos = "l",
      tl.cex = 0.6,
      tl.col = "black",
      tl.offset = 0.25,
      cl.length = 9,
      cl.cex = 0.6,
      cl.ratio = 0.25,
      family = "mono"
  )

Correlação

A maioria das variáveis independentes não apresentam uma alta correlação. No entanto, os níveis de insulina tem uma alta correlação com a glicose, o IMC tem uma alta correlação com a espessura da dobra cutânea do tríceps e a idade tem uma alta correlação com a quantidade de gravidez.

Métricas utilizadas para avalização dos modelos

Métricas utilizadas para avalização dos modelos

  • ROC AUC

  • Acurácia

  • Brier Score: \(BS = \frac{1}{N}\sum_{t=1}^N\left(f_t - o_t\right)^2\). \(N\) é o tamanho da amostra classificada, \(f_t\) é a probabilidade predita e \(o_t\) é o valor observado.

Preparação dos dados

Preparação dos dados

data_split <- initial_split(data, prop = 0.8, strata = Outcome)
train_data <- training(data_split)
test_data <- testing(data_split)

rec <- recipe(Outcome ~ ., data = train_data) |>
  step_impute_median(all_numeric_predictors()) |>
  step_mutate(across(all_numeric_predictors(), log1p)) |>
  step_mutate(
    N0 = BMI * SkinThickness,
    N8 = Age / Insulin,
    N13 = Glucose / DiabetesPedigreeFunction,
  ) |>
  step_nzv(all_numeric_predictors()) |>
  step_normalize(all_numeric_predictors())

Preparação dos dados

  • Os dados foram divididos entre teste e treino, com 20% para o conjunto de teste e estratificado pela variável dependente.

  • Os valores ausentes foram substituídos pela mediana das variáveis.

  • Foram criadas novas variáveis, idade por nível de insulina, glicose por predisposição de ter diabetes e o produto entre IMC e a espessura da dobra cutânea do tríceps.

  • Todas as variáveis foram transformadas com \(\log \left(1 + x\right)\) e normalizadas. Além disso, foi utilizado o utilizado o step_nzv, que serve para remover variáveis com variância próximo de zero.

Algoritmos e validação cruzada

set.seed(123)

cv_folds <- vfold_cv(train_data,
                     v = 10,
                     strata = Outcome
                     )

knn_spec <- nearest_neighbor(neighbors = tune(), weight_func = tune()) |>
  set_engine("kknn") |>
  set_mode("classification")

nbayes_spec <- naive_Bayes(smoothness = tune(), Laplace = tune()) |>
  set_engine("naivebayes") |>
  set_mode("classification")

lda_spec <- discrim_linear() |>
  set_engine("MASS") |>
  set_mode("classification")

qda_spec <- discrim_quad() |>
  set_engine("MASS") |>
  set_mode("classification")

logistic_spec <- logistic_reg(
  engine = "glmnet",
  penalty = tune(),
  mixture = tune()
)

dt_spec <- decision_tree(
  cost_complexity = tune(),
  min_n = tune(),
  tree_depth = tune()) |>
  set_engine(engine = "rpart") |>
  set_mode("classification")

rf_spec <- rand_forest(mtry = tune(),
                       min_n = tune(),
                       trees = tune()) |>
  set_engine(engine = "ranger") |>
  set_mode("classification")

bt_spec <- bag_tree(
        cost_complexity = tune(),
        min_n = tune(),
        tree_depth = tune()
    ) |>
    set_engine("rpart") |>
    set_mode("classification")

rf_spec <- rand_forest(
        mtry = tune(),
        min_n = tune(),
        trees = tune()
    ) |>
    set_engine(engine = "ranger") |>
    set_mode("classification")

xgb_spec <- boost_tree(
    tree_depth = tune(),
    learn_rate = tune(),
    loss_reduction = tune(),
    min_n = tune(),
    sample_size = tune(),
    trees = tune(),
    mtry = tune()
    ) |>
    set_engine(engine = "xgboost") |>
    set_mode("classification")

lsvm_spec <- svm_linear(cost = tune(),
                        margin = tune()) |>
  set_engine(engine = "kernlab") |>
  set_mode("classification")

rsvm_spec <- svm_rbf(cost = tune(),
                     rbf_sigma = tune(),
                     margin = tune()) |>
  set_engine(engine = "kernlab") |>
  set_mode("classification")

psvm_spec <- svm_poly(cost = tune(),
                      degree = tune(),
                      scale_factor = tune(),
                      margin = tune()) |>
  set_engine("kernlab")  |>
  set_mode("classification")

wf = workflow_set(
  preproc = list(rec),
  models = list(
        knn_fit = knn_spec,
    nbayes_fit = nbayes_spec,
    lda_fit = lda_spec,
    qda_fit = qda_spec,
    linear_fit = logistic_spec,
        dt_fit = dt_spec,
        bt_fit = bt_spec,
        rf_fit = rf_spec,
        xgb_fit = xgb_spec,
    lsvm_fit = lsvm_spec,
    rsvm_fit = rsvm_spec,
    psvm_fit = psvm_spec
    )
  )  |>
  mutate(wflow_id = gsub("(recipe_)", "", wflow_id))

grid_ctrl = control_grid(
  save_pred = TRUE,
  parallel_over = "resamples",
  save_workflow = TRUE
)

grid_results = wf  |>
  workflow_map(
    seed = 123,
    resamples = cv_folds,
    grid = 50,
    control = grid_ctrl
  )
autoplot(grid_results, metric = "roc_auc") +
  theme_bw() +
  labs(y = "Métrica", x = "")
autoplot(grid_results, select_best = TRUE) +
  theme_bw() +
  labs(y = "Métrica", x = "")

Algoritmos e validação cruzada

  • A validação cruzada foi realizada com 10 folds, estratificado pela variável dependente e foi criado um grid de 50 combinações entre os hiperparâmetros. Assim, foram otimizados os algoritmos a quantidade de vizinhos e o tipo de distância utilizada. Foram otimizados, para o naive bayes, o parâmetro de suavidade para o limite da classe e a correção da suavidade chamado de Laplace. Foram também otimizados os hiperparâmetros de regularização da regressão logística.

Resultado geral da otimização

Resultado na base de teste

results_acc = workflowsets::rank_results(
  grid_results,
  select_best = TRUE,
  rank_metric = "roc_auc"
  ) |>
  filter(.metric == "roc_auc") |>
  dplyr::select(wflow_id, mean, std_err, model, rank)

model_ids <- c(
  "linear_fit", "knn_fit",
  "nbayes_fit", "lda_fit",
  "qda_fit", "dt_fit",
  "bt_fit", "rf_fit", "xgb_fit",
  "rsvm_fit", "lsvm_fit", "psvm_fit")

best_sets <- map(model_ids, ~ grid_results %>%
                   extract_workflow_set_result(.x) %>%
                   select_best(metric = "roc_auc"))

names(best_sets) <- model_ids

test_results <- function(rc_rslts, fit_obj, par_set, split_obj) {
  rc_rslts %>%
    extract_workflow(fit_obj) %>%
    finalize_workflow(par_set) %>%
    last_fit(
      split = split_obj,
      metrics = metric_set(accuracy, roc_auc, brier_class
      )
    )
}

test_results_list <- map2(model_ids, best_sets,
                          ~ test_results(grid_results, .x, .y, data_split))

metrics_table <- map_dfr(test_results_list, collect_metrics, .id = "model_id") |>
  dplyr::select(model_id, .metric, .estimate) |>
  pivot_wider(names_from = .metric, values_from = .estimate) |>
  mutate(across(where(is.numeric), round, 4)) |>
  mutate(
    Modelo = c("Logistic", "K-NN",
    "Naive Bayes", "LDA",
    "QDA", "Árvore de Decisão", "Boost tree",
    "Floresta aleatória", "XGBoost",
    "RSVM", "LSVM", "PSVM"
    )
  ) |>
  relocate(Modelo) |>
  dplyr::select(-model_id)

colnames(metrics_table) <- c("Modelo", "Accuracy", "ROC AUC", "Brier Score")

metrics_table <- as_tibble(metrics_table)
metrics_table |>
  arrange(desc(`ROC AUC`)) |>
  knitr::kable()
Modelo Accuracy ROC AUC Brier Score
Boost tree 0.7857 0.8787 0.1357
LSVM 0.7922 0.8681 0.1465
LDA 0.7857 0.8652 0.1443
Logistic 0.7792 0.8635 0.1484
XGBoost 0.7662 0.8593 0.1695
PSVM 0.7792 0.8569 0.1516
RSVM 0.7792 0.8509 0.1517
Floresta aleatória 0.7662 0.8431 0.1535
K-NN 0.7792 0.8355 0.1589
QDA 0.7727 0.8238 0.1663
Árvore de Decisão 0.7922 0.8197 0.1594
Naive Bayes 0.7013 0.8128 0.1842

Resultado na base de teste

O modelo que obteve o melhor resultado foi o discriminante linear, com base na curva ROC e o brier. Na base de teste, ele obteve um ROC AUC de 0,8652, a métrica de brier 0,1443 e obteve uma acurácia de 78,57%. O que pior perfomou foi o Naive Bayes, obtendo as piores estatísticas para todas as métricas consideradas.

Hiperparâmetros selecionados dos modelos tunados

# Hiperparâmetros para o K-NN

grid_results |>
  extract_workflow_set_result("knn_fit") |>
  select_by_one_std_err(desc(neighbors), metric = "roc_auc") |>
  knitr::kable()
neighbors weight_func .config
15 gaussian Preprocessor1_Model22
# Hiperparâmetros do naive bayes

grid_results |>
  extract_workflow_set_result("nbayes_fit") |>
  select_by_one_std_err(desc(Laplace), metric = "roc_auc") |>
  knitr::kable()
smoothness Laplace .config
1.132653 3 Preprocessor1_Model32
# Hiperparâmetros da regressão logística com penalização

grid_results |>
  extract_workflow_set_result("linear_fit") |>
  select_by_one_std_err(desc(mixture), metric = "roc_auc") |>
  knitr::kable()
penalty mixture .config
0.0002121 1 Preprocessor1_Model50
# Hiperparâmetros para bagged trees

grid_results |>
  extract_workflow_set_result("bt_fit") |>
  select_by_one_std_err(desc(tree_depth), metric = "roc_auc") |>
  knitr::kable()
cost_complexity tree_depth min_n .config
1.39e-05 15 28 Preprocessor1_Model29
# Hiperparâmetros da random forest

grid_results |>
  extract_workflow_set_result("rf_fit") |>
  select_by_one_std_err(desc(trees), metric = "roc_auc") |>
  knitr::kable()
mtry trees min_n .config
5 2000 11 Preprocessor1_Model24
# Hiperparâmetros do boosting

grid_results |>
  extract_workflow_set_result("xgb_fit") |>
  select_by_one_std_err(desc(tree_depth), metric = "roc_auc") |>
  knitr::kable()
mtry trees min_n tree_depth learn_rate loss_reduction sample_size .config
5 449 8 15 0.0014225 4.9e-06 0.6510204 Preprocessor1_Model22
# Hiperparâmetros da árvore de decisão

grid_results |>
  extract_workflow_set_result("dt_fit") |>
  select_by_one_std_err(desc(tree_depth), metric = "roc_auc") |>
  knitr::kable()
cost_complexity tree_depth min_n .config
1.39e-05 15 28 Preprocessor1_Model29

Referência

Pima Indians Diabetes. Baixado através do Kaggle em: https://www.kaggle.com/datasets/uciml/pima-indians-diabetes-database.